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ABSTRACT 

Wc present initial results from "Via Lactea" , the highest resolution simulation to date of Galactic 
CDM substructure. It follows the formation of a Milky Way-size halo with Mhaio = 1-8 x 10^^ Mq 
in a WMAP 3-year cosmology, using 234 million particles. Over 10,000 subhalos can be identified 
at z=0: Their cumulative mass function is well-fit by A^(> Mgub) = 0.0064 (Afgub/Afhaio)"^ down 
to Msub = 4 X 10^ M0. The total mass fraction in subhalos is 5.3%, while the fraction of surface 
mass density in substructure within a projected distance of 10 kpc from the halo center is 0.3%. 
Because of the significant contribution from the smallest resolved subhalos, these fractions have not 
converged yet. Sub-substructure is apparent in all the larger satellites, and a few dark matter lumps 
are resolved even in the solar vicinity. The number of dark satellites with peak circular velocities above 
10 kms~^ (5 kms~^) is 124 (812): of these, 5 (26) are found within O.lrvir, a region that appeared 
practically smooth in previous simulations. The neutralino self-annihilation 7-ray emission from dark 
matter clumps is approximately constant per subhalo mass decade. Therefore, while in our run the 
contribution of substructure to the 7-ray luminosity of the Galactic halo amounts to only 40% of the 
total spherically-averaged smooth signal, we expect this fraction to grow significantly as resolution is 
increased further. An all-sky map of the expected annihilation 7-ray flux reaching a fiducial observer 
at 8 kpc from the Galactic center shows that at the current resolution a small number of subhalos start 
to be bright enough to be visible against the background from the smooth density field surrounding 
the observer. 

Subject headings: cosmology: theory - dark matter - galaxies: dwarfs - galaxies: formation - galaxies: 
halos - methods: numerical 



1. INTRODUCTION 

In structure formation scenarios dominated by cold 
dark matter (CDM), the halos of galaxies and galaxy 
clusters are assembled via the hierarchical merging and 
accretion of smaller progenitors. This process causes 
structures to virialize to a new equilibrium by redis- 
tributing energy among the coUisionlcss mass compo- 
nents. Early low resolution numerical simulations and 
simple analytical models found that the end products 
of this "bottom-up" scenario were smooth, triaxial ha- 
los. In recent years, higher resolution cosmological N- 
body simulations have modified this picture: the merg- 
ing of progenitors is not always complete, and the cores 
of accreted halos often survive as gravitationally bound 
subhalos orbiting within a larger host system. CDM 
halos are not smooth, they have a wealth of substruc- 
ture on all resolved m ass scales (e.g. Moore et al."199^ 
Klvpin et al.„ 1999: , Moore et al.ir20 01: Stochr et al. 2002 
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The amount and spatial distribution of subhalos 
around their host provide unique information and clues 
on the galaxy assembly process and the nature of the 
dark matter. The mismatch between the dozen or so 
dwarf satellite galaxies observed around the Milky Way 
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and the predicted large n umber of CDM subhalos of com- 
parab le circular velocity (jMoore et al.lll999l : lKlvpin et al.l 
|1999( ) has been the sulDject of many recent studies. So- 
lutions involving feedback mechanisms that make star 
formation in smal l halos very inefficient offer a possi- 
ble wav out (e . g. iBuUock et all [2OOOI: ISomervilld I2OO2I: 



iKravtsov etall 120041: iMoore et al.ll2006f) . Other models 
have attempted to solve the apparent small-scale prob- 
lems of CDM at a more funda mental level, i.e. by re- 
ducin g small-scale power (e.g. iKamionkowski fc Liddld 
|2000( ). Even if most dark matter satellites have no 
optically luminous counterparts, the substructure pop- 
ulation may be detectable via flux ratio anomalies in 
strong gravitational lenses (e.g. iMetcalf fc Madaull2001l : 
[Chiba 2003). through i ts effects on stel lar streams 
(e.g. Ilbata et all 120021 : iMaver et al] l2002l ). or possi- 
bly vi a 7-rays from d ark matt er annihilation in their 
cores (jBergstrom et al.l [T999: Ca lcaneo-Roldan fc Moorel 
20001: IStoehr et all l2003t iColafrancesco et all l2006l : 
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The possibility of observing the fingerprints of the 
small-scale structure of CDM hinges on the ability of 
subhalos to survive the hierarchical clustering process as 
substructure within the host. This is in turn particularly 
sensitive to resolution issues, as subhalos with numeri- 
cally softened central densities arc more easily disrupted 
by tidal forces. In this paper we present initial results on 
halo substructure from a new N-body cosmological sim- 
ulation of unprecedented dynamic range: dubbed "Via 
Lactea" , it resolves the galaxy forming region of a Milky 
Way-size halo at z = with over 200 million particles. 
This is an order of magnitude more than used in previ- 
ous simulations. The run was completed in 320,000 CPU 
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Zj Zf e Wp mp 

48.4 90 pc 2.34 x 10» 2.09 x IC* Mq 
48.4 0.0 378 pc 6.50x10^ 5.64 x 10^ M© 

Note. — Initial and final redshifts Zi and Zf, (spline) softening 
length e, total number A^p and mass mp of dark matter particles 
for the Via Lactea simulation (top) and a lower resolution run 
{bottom). Force softening lengths t are constant in physical units 
back to 2: = 9 and constant in eomoving units before. 

hours on NASA's Project Columbia supercomputer, cur- 
rently one of the fastest machines available. 

2. SIMULATIONS 

The simulation was performed with PK DGRAV 
(lStade]|[200lt iWadslev. Stadel fc Quin^[200l and em- 
ployed multiple mass particle grid in itial conditions gen - 
crated with the GRAFICS package ()Bertschingedl2001f ). 
The high resolution region was sampled with 234 mil- 
lion particles particles of mass 2.1 x 10'' Mq and evolved 
with a force resolution of 90 pc. It was embedded 
within a periodic box of eomoving size 90 Mpc, which 
was sampled at lower resolution to account for the large 
scale tidal forces. We adopted the best-fit cosmolog- 
ical p arameters from the WMAP three-year data re- 
lease (jSpergel et al.ll2006[) : = 0.238, f^A = 0.762, 
Ho = 73 kms-i Mpc-\ n = 0.951, and erg = 0.74. The 
lower values of a% and n compared to 1-year WMAP 
results have the effect of delaying structure formation 
and reducing small-scale power. The simulation was cen- 
tered on an isolated halo that had no major merger after 
z = 1.7, making it a suitable host for a Milky Way- like 
disk galaxy. In this work we focus on results obtained 
at z=0, the formation history of the Via Lactea halo and 
its substructures will be presented in a subsequent paper. 
The host halo mass at z = is Mhaio = 1-77 x 10^^ Mq 
within a radius of rvir = 389 kpc (We define " rvir" as the 
radius within which the enclosed average density is 200 
times the mean matter value). We used adaptive time- 
steps as short as (13.7/2x10^) Gyr = 68, 500 yr satisfying 
the condition Ai < 0.2-y/e/a, where e is the force soften- 
ing length and a is the local acceleration. This criterion 
is insufficient to ensure convergence at small radii, where 
the dynamical time becomes too short and runs using the 
\ftja criterion start to underestimate matter densities. 
In comparable runs this problem was observed to to give 
10% too lo w local densities at 6 x 10^pcrit(2^ = 0), see 
Figure 1 of iDiemand et al1l2005l . This density is reached 
at 0.0034 Tvir in the Via Lactea run, therefore we set 
'"conv = 0.0034 Tvir. The local densities (but not neces- 
sarily the logarithmic slopes) should be correct to within 
10% beyond 0.0034 rvir, but detailed convergence tests 
using galaxy scale halos would be needed to confirm this. 
To resolve the density profile further in than 0.0034 rvir, 
adaptive time-steps that sc ale with dynamic al timescale 
should be adopted instead (jZemp et al.ll2006[ ) . The force 
calculations in the tree algorithm PKDGRAV includes 
multipole moments up to hexa-decapole to reach high 
force accuracy. To check for numerical convergence of 
substructure properties we resimulated the same halo 
(plus a few neighboring galaxy halos) with 27 times more 
massive particles. The parameters of the two simulations 
are summarized in Table [21 

Figure [1] (upper panel) shows the logarithmic slope of 
the spherically-averaged density profile, dlogp/dlogr, 



for our Via Lactea halo, plo tted versus radius. The 
fitting formula proposed by iNavarro. Frenk. fc White! 
()1997() with scale radius of 24.6 kpc provides a reasonable 
approximation to the density profile down to our con- 
vergence radius rconv Within the region of convergence, 
deviations from the best-fit NFW matter density are typ- 
ically less than 10%. From 0.025 rvir down to rconv Via 
Lactea is actually denser than predicted by the NFW for- 
mula. Near rconv the density approaches the NFW value 
again while the logarithmic slope is shallower (—1.0 at 
?'conv) than that the NFW fit. While we have checked 
that these conclusions are robust to variations in the 
binning used to construct the density profile, they might 
be affected by numerical flattening near our convergence 
radius. We pl an to run n ew simulations with improved 
time-stepping (jZemp et al..,2006i ) to address this issue. 

3. SUBSTRUCTURE IN GALAXY HALOS 

The wealth of substructure resolved in our Via Lactea 
run is clearly see in Figure [2l About 10,000 surviving 
satellites, some sub-subhalos (a few examples are shown 
in Fig. [3]), and even a couple of dark matter clumps at 
less than 8 kpc from the Galactic Center are now re- 
solved. 

3.1. Subhalo identification 

We have constructed substructure catalogs using 
the phase -space friends-o f -friend s algorithm 6DF0F de- 
scribed in lDiemand et al.l (|2006D . We adopt a space fink- 
ing length of 0.28 physical kpc (0.033 of the mean particle 
separation at z = 0), a velocity space linking length of 16 
physical kms^^ (0.1 of the mean ID velocity dispersion 
of the host halo), and keep all groups with at least 16 
members. Around their centers, spherical density pro- 
files (using 50 logarithmic bins) were constructed: the 
bins reach out to rvir in the case of isolated halos or out 
to the center of the nearest larger halo in the case of 
subhalos, i.e. all structures are allowed to contain their 
own, smaller substructures. The resulting circular veloc- 
ity profiles were then fitted with the sum of an NFW pro- 
file and a constant density background. The radial range 
used in these fits starts at 3 force softening lengths, i. e. 
0.27 kpc in the Via Lactea run, and ends either at the 
largest radial bin stored during at runtime or further in if 
the circular velocity grows linearly (within 10%) over five 
consecutive bins. Figure 2] shows nine example subhalo 
circular velocity profile fits. A subhalo tidal radius rt 
was defined as Psuhift) = 2 x pbg, where pbg is the local 
matter density of the host halo. This is the tidal radius 
of an isothermal (psub cx r~^) satellite on a circular orbit 
within an isothermal host. We define the subhalo mass 
to be all the mass within rt, i.e. we did not perform an 
unbinding procedure. The speed of 6DF0F and its par- 
allel implementation into PKDGRAV allowed us to use 
it at run-time with negligible computational cost (com- 
pared to the gravity calculations) on all 200 snapshots 
that were stored (separated by 68.5 Myr). 

The 6DF0F subhalo list was checked in various regions 
of the Via La ctea halo against a list obtained with SKID 
(|Stadelll200ll ). The results compare very well down to 
about 40 particles per halo except for the very largest 
subhalos, where the SKID algorithm tends to underesti- 
mate the masses and is unable to properly account for 
the presence of resolved sub-substructure. Detailed evo- 



Dicmand, Kuhlen, & Madau 



3 




Fig. 1. — Top: logarithmic slope of the density profile of our Via Lactea run, as a function of radius. Densities were computed in 50 radial 
logarithmic bins, and the local slope was determined by a finite difference approximation using one neighboring bin on each side. The thin 
line shows the slope of the best-fit NFW profile with scale radius of 24.6 kpc. The vertical dotted line indicates the estimated convergence 
radius: local densities (but not necessarily the logarithmic slopes) should be correct to within 10% outside of this radius. Bottom: the 
residuals in percent between the density profile and the best-fit NFW profile, as a function of radius. 



lutionary tracks of subhalo properties and density profiles 
will be published in a subsequent paper. 

3.2. Subhalo abundance 

In Figure [5] we present the cumulative mass func- 
tion for our entire subhalo populationw ithin rvir as 
well as for the subpopulation within the inner 0.1 rvir- 
The former is well- fit by a power-law A^(> Mguh) = 
0.0064(A/™b/Mhaio)~' above Af^b = 4 x lO^ Mq, as 
shown by the solid line (formally, the best-fit slope 



in the mass range Msub = 4 x 10'^ - 4 x 10^ M© is 
—0.97 ± 0.03). This slope implies an equal mass per 
subhalo mass decade, making the total mass in substruc- 
ture quite dependent on resolution. Note how numeri- 
cal resolution effects start to flatten the distribution al- 
ready below masses corresponding to about 200 particles 
(i.e. below M^uh = 4 x 10'^ M© and M^uh = 10® M© 
for the Via Lactea and the lower resolution run, respec- 
tively). The fraction of the host halo mass in subha- 
los is fsuh = 5.3% (0.87%) within rvir (0.1 7Vir). The 




Fig. 2. — Projected dark matter density-squared map of our simulated Milky Way-size halo ("Via Lactea") at the present epoch. The 
imago covers an area of 800 X 600 kpc, and the projection goes through a 600 kpc-deep cuboid containing a total of 110 million particles. 
The logarithmic color scale covers 20 decades in density-square. 
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Fig. 3. — Projected dark matter density-square map of the four most massive subhalos within the simulated Milky Way host at the present 
epoch. Sub-substructure is clearly visible. Only dark matter particles within the tidal radius rt are used for the projections. Clockwise 
from top left: (Msuh,rt,rv^^J = (9.8 X 10^ Mq, 40.1 kpc, 7.6kpc), (3.7 X 10^ Mq, 33.4kpc, 4.0kpc), (3.0 X lO^ , 28.0 kpc, 4.9 kpc), and 
(2.4 X 10^ Mq, 14.7kpc, 6.1 kpc). The mean subhalo densities within the tidal radius (in units of the cosmic background dark matter 
density) are 1002, 654, 904, and 4950, respectively. These values are related to the local matter density of the host (72, 46, 59 and 397 in 
the same units), and correlate only weakly with the subhalo distance from the Galactic center (345, 374, 280 and 185 kpc). 
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Fig. 4. — Examples of subhalo circular velocity profiles {crosses). The fitting functions {dashed lines) are combinations of an NFW 
subhalo {dotted lines) and a linear contribution form a constant background density {solid lines). Stated above each panel are the subhalos 
distance form the galaxy center, tidal radius, tidal mass and the local background density. 



contribution from subhalos below 10® Mq is significant: 
/sub(Msub < 10® Mo) = 2.3% (0.15%). The total sub- 
structure mass fractions in the lower resolution run are 
substantially smaller, /sub = 2.6% (0.11%). This is in 
part because the low resolution realization happens to 
contain slightly fewer massive satellites at z = 0, and 
in part because of the reduced number of resolved small 
subhalos. 

Figure [6] shows the cumulative peak circular velocity 
(Knax = max{-y/G'Msub(< function for the en- 

tire subhalo population within Tvir as well as for the in- 
ner sub-population. The former is well-fit by a power- 
law N{> V^nax) = (l/48)(K„ax)/Knax,host)"^ abovC 

5 kms~^, as shown by the dotted line. This is in good 
agreement with the mean a bundance obtaine d from a 
large sample of host halos bv lReed et al.l (|2005( l after ac- 
counting for the effect of the different cosmology.^ Note 
how numerical effects start to flatten the velocity func- 
tions below 5 kms^^ and already below 15 kms~^ for 
the lower resolution run. 

3.3. A "missing inner satellites" problem? 

^ IReed et al.l 1)20051 ') simulations were normalized using as = 1-0. 
Decreasing trg to 0.74 (the value used here) lowers the subhalo 
veloc ity function by ab out a factor two (see IZentner &: Bulloctd 
[2003 and the erratum to IDiemand eTaIll2004lJ) . 



The number of dark matter subhalos with peak circu- 
lar velocity larger than 10 kms~^ is 124 (compared to 
74 in the lower resolution run) . This is an order of mag- 
nitude larger than the dozen or so dwarf satellite galax- 
ies of comp arable circular velocity observed arou nd the 
Milky Way (jMoore et al.lll999l : lKlvpin et al.lll999f ). a well 
known mismatch often referred to as the "missing satel- 
lites problem" . The WMAP 3-year cosmological param- 
eters do not alleviate the discrepancy significantly, and 
with our very high resolution Via Lactea the difference 
between the predicted number of CDM subhalos and the 
satellites actually observed has become even more pro- 
nounced. Indeed, we find a new "missing inner satellites 
problem" within 0.1 rvir = 39 kpc, a region that appeared 
practically smooth in previous lower resolution simula- 
tions. In this inner region we now resolve 5 (26) dark 
satellites with Knax > 10 kms~^ (Knax > 5 kms~^), 
compared to only one known Milky Way dwarf galaxy, 
Sagittarius at a distance of 24±2 kpc (|Mateolll998f l. Even 
when the suspected stripped cores of dwarf spheroi dals 
from the recent compilation of lvan den Berghl (j2006[ ) are 
included, the numbers of observed inner satellites remain 
well below the predicted number of substructures in the 
inner halo. 

Recent Local G roup models (jKravtsov et al.l l2004l : 
iMoore et all l2006[ ) were based on simulations like the 
lower resolution run used in this work. The overall abun- 
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dance of subhalos does increase with resolution in the 
relevant V,nax range, but not dramatically, i.e. the to- 
tal abundance of satellite galaxies in these models seems 
not to be very affected by their lower numerical resolu- 
tion. But what would these two models predict about 
satellite gal axies within the newl y resolved inner sub- 
halos? The iKravtsov et al.l (|2004[ ) model allows galaxy 
formation only in subhalos that were relatively massive 
(> 10^ Mq) before the y fell into the gravit ational poten- 
tial of the host. In the lMoore et al.l (|2OO60 scenario only 
the earliest forming halos above the atomic cooling mass 
at z > 12 become proto-galactic building blocks (many 
of these systems remain well below 10^ M0 at all times). 
By z = they have built up a halo of stars, globular clus- 
ters and a few surviving dwarfs. We followed the assem- 
bly history of the inner subhalos backwards in time and 
found that only two of the m lie above the (time depen- 
dent) minimum mass from IKravtsov et all l|2004D . The 
same two systems are also the only o nes that forni early 
enough to become luminous in the |Moore et al.l 120061 ) 
scenario. In the inner O.lrvir both models thus predict 
two dwarf galaxies, in good agreement with the numbers 
observed around the Milky Way and M31, and leave the 
other three inner subhalos above 10 kms~^ dark. There- 
fore the new, missing inner satellites problem seems to 
be resolvable in the same way as the well known Local 
Group wide problem. The overabundance of inner sub- 
halos might not be a problem for CDM, but it leads to 
the new and interesting prediction of having a number 
of relatively large, dark CDM subhalos orbiting in the 
inner halo, i.e. in the same region where the Milky Way 
galaxy is located. 



3.4. Mass fractions in substructure and gravitational 

lensing 

Standard smooth gravitational lens models have 

difficulties in explaining the relative flu xes of 

mult iply- imaged quas ars (|Mao fc Schneider! Il998l : 
iMetcalf fc Zhao I l2002t l. The discrepancy between 
the predicted and observed flux ratios is commonly 
referred to as the "anomalous flux ratio problem" , and 
dark matter substructure within the lens halo is one 
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fraction, lensing observations are sensitive to the mass 
fraction in substructure projected through a cylinder of 
radius 5-10 kpc around the lens center. Figure [7] depicts 
the fraction of the host halo mass within a sphere of 
radius r that is bound up in substructure, /sub(< r), 
as well as the substructure mass fraction in cylindrical 
projection of radius R, /sub(< R), measured in our Via 
Lactca simulation. The radial distribut ion follows the 
subhalo number density profile given in iDiemand et al.l 
(|2004b( l. i.e. it is more extended than the overall mass 
distribution. In the vicinity of i? = 10 kpc, the surface 
mass density for all subhalos with Mgub < 10^ M0 can 
be approximated as /sub(< R) = 0.002 (i?/10kpc). 

Our total projected surface densities are consistent 
with, but on the low side of, estima tes from semiana- 
lytic models (jZentner fc Bulloc3l2003l ) , although we find 
significantly larger contributions from the smallest sub- 
halos at large projected radii. The total projected sur- 
face densities is lower than the few percent value that 
seems to be required to explain the anomalous flux ratios 
(jDalal fc Kochanek|[200l iMetcalf et an[200l . Whether 
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Fig. 6. — Cumulative peak circular velocity function for all subhalos within r^ir {upper panel) and for the subpopulation within the 
inner O.lrvir {lower panel). Solid lines without error bars: Via Lactea run. Dashed lines: lower resolution run. Solid line: power-law fit 
N{> Vmax ) = (1/ 48) X ( Vmax/Vmax, host) ^- Solid lines with error bars: observed number of dwarf galaxy satellites around the Milky 
Way (as in lMateol l[T§9l ) : IKlvpin et al.l l[T99l ). plus the recently discovered Ursa Major dwarf HWillman et al. Il2005l : IKlevna et al.|[2005^ '). 
within r^ir and within O.lrvir (only Sagittarius). The vertical error bars show the Poissonian y/N scatter. 
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r [kpc] , R [kpc] 

Fig. 7. — Cumulative mass fraction of substructure, shown spherically-averaged as a function of radius r {dashed lines) and in cylindrical 
projection as a function of projected radius R [solid lines). Averages over 1000 random projections arc shown, and only subtructures within 
Tvir are included. Thick lines: all subhalos. Thin lines: subhalos with M^^^, < (10^ Mq, 10* Mq, 10® Mq) from bottom to top. 
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this constitutes a serious disagreement with CDM predic- 
tions is unclear. Here we just notice that our simulation 
likely suffers the numerical "over-merging" problem be- 
low 4 X 10^ M©: if the mass spectrum of substructure 
we measure, dN/dMsuh oc M^, extended all the way 
down to lO'' M0, for example, then the mass fraction 
in subhalos with iC M© < Msub < 4 x 10^ M© would 
be comparable to that in all higher mass satellites with 
4 X 10^ < Msub < 10^ Mq identified in our run. We also 
have to expect some "over-merging" in the dense, inner 
halo: Via Lactea is able to resolve a few subhalos within 
10 kpc of the galactic center (3 systems at z = 0), but 
the true subhalo abundance might be significantly higher 
in this region. 

4. DARK MATTER ANNIHILATION SIGNAL 

Our simulation shows that, in the standard CDM 
paradigm, galactic halos should be filled with at least 
tens of thousands subhalos that have no optically lumi- 
nous counterpart. If the dark matter is in the form of a 
supersymmetric particle produced in the early Universe 
like the neutralino, then substructure may be lit up by 
the annihilation of such particles into 7-rays. Since the 
annihilation rate is proportional to density squared, the 
predicted flux depends sensitively on the dumpiness of 
the mass distribution. ^ Within the Milky Way, nearby 
subhalos may be among the brightest sources of annihi- 
lation radiation, and could be detectable by the forth- 
coming Gamma-Ray Large Area Telescope (GLAST). 

Annihilation luminosities based on densities measured 
directly in coUisionless N-body simulations ^ suggest 
that bound substructure may boost the signal from 
individual halos by only a small f actor compared to 
a smooth spherical density profile (jStoehr et al.l l2003l : 
iDiemand et al.l l2006( l , and that the cnhancenement is 
dominated by the most massive subhalos. The measured 
contribution of resolved substructure must be affected 
by numerical resolution, however, since higher resolu- 
tion simulations would b e able to resolve higher den- 
sities in subhalo centers (jKazantzidis et ani2004[ ). In- 
deed, analytic calc ulations tend to find larger subs truc- 
ture boost factors (jCalcaneo-Roldan fc Mooril2000f ) and 
a signal that is dominate d by small subhalos instead 
(|Colafrancesco et al.l I2006D . In Figure M we show the 
mass dependence of the annihilation luminosity for indi- 
vidual subhalos in our two Milky Way simulations. The 
signal from the i'^ subhalo is proportional to 

5z = / PsuhdV, V ftmp, (1) 

where pj is the density of the j*'^ particle, and {Pi} is 
the set of all particles belonging to halo i. In the Via 
Lactea run we find ^sub oc Afgub, i-e. a signal-to- mass 
ratio that is approximately independent of mass. Given 

^ Annihilation is too slow to alter the mass distribution on the 
scales simulated here. For an NFW subhalo it would erase the cusp 
only on scales smaller than a micro-parsec and this would not lower 
the total subhalo luminosity significantly. 

In this work we neglect the still somewhat unclear effects of 
galaxy formation on the inner dark matter density. Substantial 
gains are possible, e.g . the contracted NFW halo "NFWc" in 
IMambrini et all I I2006I) . rescaled to match Via Lactea, gives an 
eight times larger host halo signal than the corresponding NFW 
profile before contraction. 



our substructure abundance of dN/dMsuh k M^^, this 
implies an annihilation luminosity that is approximately 
constant per decade of substructure mass, as the Figure 
shows. A comparison with our lower resolution simula- 
tion indicates that more than 500 dark matter particle 
per subhalo are needed to produce a signal that is of the 
right order of magnitude. Below this scale (i.e. below 
about 10^ Mq in Via Lactea and 3 x 10® Mq for the lower 
resolution run) the luminosity decreases quickly due to 
insufficient resolution. 

To test the robustness of the above results we have 
tried two alternative ways to compute matter densities 
and annihilation luminosites besides the standard SPH- 
kernel. The first is based on densities measured in spheri- 
cal bins around each subhalo, while the second uses NFW 
fits to these density profiles. All three methods agree 
with an approximately constant total signal per decade 
of substructure mass above a mass scale of 500 parti- 
cles per subhalo. However, the spherical density profile 
estimate falls below the NFW estimate even for mas- 
sive subhalos and the difference increases with decreas- 
ing subhalo mass . This is consistent with the result of 
iKazantzidis et aP (|2004[ ) that well resolved subhalos re- 
main as cuspy as halos in the field, while under-resolved 
structures have relatively large, artificial constant den- 
sity cores. Annihilation luminosities based on local SPH 
densities, on the other hand, are larger than the spheri- 
cal signal for massive subhalos because of resolved sub- 
substructures, triaxiality, and other inhomogenities that 
get averaged out in the spherical estimates. Smaller sub- 
halos, containing a few thousand particles, appear rather 
smooth and spherical and have a large artificial core. It 
is therefore not surprising that in this range the SPH es- 
timate lies above the spherical density profile estimate, 
but below the NFW signals. Below about 3000 particles 
(6 X 10 Mq) the NFW total signal estimate decreases 
and approaches the other two estimates, presumably be- 
cause it too starts to be affected by the finite resolution: 
with this low number of particles, halo density profiles 
arc artificially flattened by numerical effects out to a ra- 
dius of about N-'^/^rt = 0.07 (Dicmand ct al. 2004a). 
The scale radii of many subhalos would lie in the affected 
range and the resulting apparent scale radii of the NFW 
fits come out too large. This leads to underestimated 
scale densities and 7-ray luminosities. Above 3000 par- 
ticles, the total signal estimate based on NFW fits hints 
at a contribution to the total signal that increases with 
decreasin g subhalo mass in agreeme nt with the analytic 
model by lColafrancesco et al.l (|2006( l. 

At our highest resolution we measure a total neutralino 
7-ray luminosity from the host halo that is a factor of 2 
higher than the total spherically-averaged smooth sig- 
nal: 41% of this increase is associated with substructure, 
the other 59% is due to other deviations fr om spherical 
sym metry. At ten times lower resolution IStoehr et al.l 
(j2003) found an increase from subhalos of 25%, and in 
our 27 times lower resolution run we measure only a 
12% increase. Since the subhalo signal appears to be 
dominated by the smaller clumps, both the total and 
the subhalo signal have not converged yet and might be 
much larger for a real ACDM halo: assuming a con- 
stant boost factor of 0.17 per decade in subhalo mass 
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Fig. 8. — The annihilation signal of individual subhalos in units of the total luminosity of the spherically averaged MW host. Crosses: 
Via Lactca run. Circles: lower resolution run. The plotted range begins at a subhalo mass corresponding to 240 particles in Via Lactea and 
only 9 particles in the lower resolution run. Densities are estimated using the SPH-kernel over 32 nearest neighbors. Solid lines: sliding 
average {lower curve) and total signal over one decade in mass [upper curve for the Via Lactea run. Dashed lines: same for the lower 
resolution run. Dash-dotted lines: Alternative estimate of the total subhalo signal in the Via Lactea run based on densities measured in 
spherical bins around the subhalo centers (lower curve), and same (u ppper curve) based on NFW fits to these density profiles. Dotted line: 
mass dependence of the total signal from [Colafrancesco at aTl I I2006I) . see text for details. 



from 10^^ Mq® to 10^" Mq would give a total enhance- 
ment factor of 3, i.e. the substructure signal would be 
twice that of a smooth spherical host halo. Allowing the 
same boost for each of the subhalos (0.17 per decade from 
10~^ Mq to 0.01 Mgub) in a self-similar manner, the total 
boost for a galaxy halo is approximately 13. A contribu- 
tion th at grows with decrea s ing su bhalo mass is not ex- 
cluded: IColafrancesco et all (|2006[ ) found such a behav- 
ior using an analytic model with a dN/dMsuh oc 
subhalo mass function, subhalo concentr ations propor- 
tional to those given for field halos in iBuUock et al.l 
(|2001h . and somewhat different cosmological parameters 
{^M = 0.281, D.K = 0.719, n = 1.0, and erg = 0.89). Fig- 
ure [8] shows that the results from the high resolution run 
are consistent with their mass dependence, especially af- 
ter considering that numerical effects tend to reduce the 
substructure signal and that their importance is gradu- 
ally gro wing with decreasing subh alo mass. Integrating 
over the IColafrancesco et al.l ()2006l ) model normahzed to 

* Micro-subhalos below 10^^ Mq appear to contribute little to 
the annihilation signal of their host IjPiemand et al.ll2006l 1. 



our result at large subhalo masses gives a boost factor 
of about 8. The increase due to sub-substructures adds 
another factor of about ten, it is more dramatic than in 
the constant S/M case. We plan to constrain this factor 
further by using analytic models based on the detailed 
subhalo and sub-subhalo properties in a future paper. 

5. 7-RAYS FROM MW SUBSTRUCTURE 

The Large Area Telescope (LAT) aboard the soon-to- 
be-launched GLAST satellite may be able to detect 7- 
rays from D M annihilations origina ting at the Galactic 
center (e.g. iManibrini et al.ll2006l ). in a subhalo (e.g. 
iKoushiappas et aTl 2004") or even in a very nearby micro- 
subhalo (jKoushiapp as 2006). The LAT has an effective 
area of more than 8000 cm^ , a field of view greater than 2 
steradian, and sub-degree angular resolution at energies 
greater than 1 GeV. GLAST is expected to operate for 
more than 5 years, and will conduct an all-sky survey 
during this time. 

We have constructed all-sky maps of the expected an- 
nihilation 7-ray flux reaching a fiducial observer located 
8 kpc from the center of the Milky Way halo. Each DM 
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Fig. 9. — All-sky map of the DM annihilation flux ^/K (Mgkpc ^ sr ^) in our Via Lactea halo, for an observer located 8 kpc from 
the Galactic center. The insets show zoom-in's of a 40° X 40° region around the anticenter (left) and the brightest subhalo (right). 



particle was assigned an annihilation 7-ray flux 



where di is the distance from the observer to the i par- 
ticle. The constant of proportionality K captures the un- 
certain particle physics and is equal to N^{av) 0/(2771^), 
where Nj is the number of photons produced per annihi- 
lation, {av)o is the annihilation rate at zero temperature, 
and rriy, the restmass of the DM particle. In this work 
we consider only the dependence of the signal on the 
macroscopic DM density distribution, and present fluxes 
in units of but for any given particle physics model our 
results can be converted to physical fluxes by multiplying 
by the appropriate K. We also calculate the angular scale 
subtended by each DM particle, A6i = hi/di, where hi is 
defined to be half the radius that encloses the nearest 32 
neighbors. We then bin up each particle's contribution 
in angular bins (pixels) of size AO, typically chosen to be 
0.1°x0.1° RsSxlO^^sr. Each particle contributes to all 
pixels within its angular scale, weighted by a projected 
SPH smoothing kernel. The weights are normalized such 
that the sum over all covered pixels equals the particle's 
total flux Fi. After binning all particles, we divide each 
pixel by AO, resulting in a map of flux per solid angle <&. 
This 2D image is mapped onto a Hammer- Aitoff equal 
area projection, with the coordinate system rotated such 
that the halo center lies at the center of the projection. 

The image of Figure [5] depicts the resulting DM anni- 
hilation all-sky map. As expected the signal is strongest 



towards the center, with a flux about four orders of mag- 
nitude larger than the diffuse contribution from the anti- 
center. Several substructures are clearly visible, both 
towards the halo center and in the opposite direction 
(towards the edges of the image). The two insets show 
in greater detail a 40° x 40° region around the anti- 
center and around the brightest subhalo. Because the 
diffuse background noise is stro ngest near the cente r due 
to Galactic continuum emission (jStoehr et al.l[2003h . it is 
likely that substructure will be more readily detectable 
at anti-center, and/or at higher latitudes. 

Compared with the density squared projection in Fig- 
ure [2] however, it is striking how little substructure is ac- 
tually visible in the all-sky map. This is due to a number 
of causes. At 8 kpc (w 0.02 rvir) the observer is located 
relatively close to the halo center. The mean density at 
this distance is 6.0 x lO^pcrit, and only 4009 of our subha- 
los have central densities greater than this. Of course this 
is largely a numerical effect, since subhalos are resolved 
with much fewer particles than the host. Assuming a uni- 
versal density profile and higher subhalo concentrations, 
subhalos would be denser than their host halo at compa- 
rable r/rvir- Moreover Via Lactea still resolves only 32 
subhalos closer than 20 kpc from the observer. In our all- 
sky map distant subhalos are hard to detect because both 
their fluxes and angular sizes on the sky are reduced. In 
Figure [TO] we have plotted cumulative distribution func- 
tions of A9 and log(Fsub/-FcGntor), as well as a scatter 
plot of the two quantities against each other. Here A9 
is the projected angular size of the subhalo's tidal ra- 
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Fig. 10. — Top: Cumulative distribution functions of the ratio of subhalo flux to central flux Kjub/^contcr (solid line) and of subhalo 
angular size A8 (dotted line), for one random observer position (the same one as in Figuro|9]l. Bottom: A8 plotted against Fsub/i^contcr 
for all sublialos. Only subhalos within rvir were used in both plots. 



dius, i^sub is the total flux from all particles within each 
subhalo's tidal radius, and i^contcr = 4.7 x lO^'^ Mgkpc^^ 
is the flux from all pixels within 1° of the Galactic cen- 
ter. We have only used subhalos within the host halo's 
Tvir- The median angular size of all subhalos is 0.48°, 
and only 5 subhalos have a flux greater than 1% of the 
halo center. The brightest subhalo, at (longitude, lat- 
itude) = {^,9) « (60°, 10°), has i^sub/Fcontcr = 0.094. 

This subhalo has a mass of ~ 1.3 x 10 Mq, a peak cir- 
cular velocity of 40 kms~^ and lies at a distance of 31 



kpc from the Galactic center. It is very unlikely that 
our Milky Way halo would contain such a massive DM 
subhalo so close to the center, since it would presumably 
host a dwarf galaxy ^. The closest such object in the Lo- 
cal Group is the Small Magellanic Cl oud, it has a com- 
parab le mass and velocity dispersion ([Harris fc Zaritskvl 
|200(3[) . but lies at a distance of w 50 kpc. It would thus 
appear 1.7 times smaller and 2.8 times fainter than the 

^ This subhalo would h ost a dwarf galaxy according to both 
models discussed in Section 13.31 
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brightest subhalo in Figure [9] 

6. CONCLUSIONS 

We have reported in this paper first resuhs from the 
highest resolution simulation to date of Galactic DM sub- 
structure. Our simulation consists of over 200 million 
DM particles and follows the evolution and formation of 
a Milky Way-size halo with Mhaio = 1.8 x lO^^ 

in a 

WMAP 3-year cosmology. Here we summarize our main 
results: 

• At the present epoch we resolve approximately 
10,000 subhalos, about one order of magnitude 
more than in any previous simulation of a Milky 
Way halo. Our resolution is sufficient to even re- 
solve a few subhalos within the solar circle. In sev- 
eral of our more massive subhalos we are able to 
identify sub-substructure. 

• The cumulative subhalo mass function is consistent 
with a (A/sub/-Mhaio)~^ power law down to Mguh = 
4 X 10^ M0 (200 particles), implying an equal mass 
contribution per decade of subhalo mass. As such 
the total mass fraction in substructure is going to 
be resolution dependent, and our value of 5.3% 
is likely to be a lower limit. Extrapolating the 
dN/dMsuh oc subhalo mass function down 
to 1 Mq, for example, would lead to a substructure 
mass fraction of 20%. The fraction of surface mass 
density in resolved substructure within a projected 
distance of 10 kpc is 0.3%, again most likely an 
underestimate. 

• The cumulative subhalo velocity function is also 
well-fit by a power law: N{> Vmax) = 
(l/48)(Vmax)/Vmax,host)"^ dowu to 5 kms-i. We 
find 124 subhalos with peak circular velocities 
greater than 10 kms~^, about one order of mag- 
nitude more than the number of dwarf satellite 
galaxies of comparable circular velocity. The lower 
CTg and n in the WMAP 3-year cosmology do not 
appreciably alleviate this "missing satellite prob- 
lem" . Furthermore we report a new "inner miss- 
ing satellite problem" concerning the inner 10% 
of Tvir- Whereas previous simulations have found 
this region to be practically smooth, our simula- 
tion reveals the presence of 5 subhalos with Vmax > 
10 km s^^ and distances less than 39 kpc, compared 
to only one kn own Milky W ay satellite, Sagittarius 
at 24 ± 2 kpc (|Mateolll998l ). A prchminary analy- 
sis suggests that this number of inner satellites is 
consistent with models in which dwarf galaxy for- 
mation is limited to subhalos sufficiently massive 
([Kravtsov et al.l 120041 ) or forming in rare density 
peaks selected at z > 12 (|Moore et al.ii2006) . 



• The total DM annihilation luminosity of our sim- 
ulated halo is a factor of 2 larger than the lumi- 
nosity derived from a smooth spherically averaged 
density profile. 41% of this increase is due to re- 
solved substructure and the remaining 59% can be 
attributed to other deviations from spherical sym- 
metry. The luminosity of individual subhalos is 
linearly proportional to their mass, which together 
with the mass function implies an equal con- 

tribution to the total subhalo luminosity from each 
decade in subhalo mass. Just like the total sub- 
halo mass fraction, the total subhalo annihilation 
luminosity fraction is thus limited by numerical res- 
olution. An extrapolation down to 10~^ M0 would 
result in a factor of 3 increase in luminosity. Sig- 
nificantly larger factors due to a signal per sub- 
halo mass decade which rises slowly with decreas- 
ing mass and/or similar gains within subhalos due 
to sub-subhalos cannot be excluded. 



• We have constructed an all-sky map of the expected 
annihilation 7-ray flux reaching a flducial observer 
at 8.0 kpc. The signal is strongest in the center and 
falls off by about 4 orders of magnitude towards 
the anticenter. With our finite numerical resolu- 
tion only a small number of subhalos are bright 
enough to be visible against the background from 
the smooth density field surrounding the observer. 
We expect actual physical DM subhalos to have 
higher central densities, making them more lumi- 
nous and possibly detectable in regions sufficiently 
far from the center. The brightest subhalo in our 
simulation has a mass comparable to the SMC, is 
located 30 kpc from the center, and has a fiux equal 
to 10% of the fiux from the central 1°. 
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